function stats = StockRetStats(parms,gesol,prob_Nxz,Phi_Nxz,stats_in,FLAG_RUN_PARALLEL)
% note: gesol.stock.price_cons_unlev is the cum-dividend price

    if nargin<6; FLAG_RUN_PARALLEL = false; end

    stats = stats_in;
    
    NN = parms.NN; Nx = parms.Nx; Nz = parms.Nz;

    price_ex_div = gesol.stock.price_cons_unlev - gesol.C;
    price_cum_div = gesol.stock.price_cons_unlev;

    % 1 month returns
    % note: multiply by 12 to annualize, multiply by 100 to convert to
    %       percentage terms.
    exp_ret_stock_t_1mo = reshape(Phi_Nxz*price_cum_div(:),[NN Nx Nz])./price_ex_div;
    rf_t_1mo = 1./gesol.real_bonds.prices(:,:,:,1);
    cond_vol_ret_t_1mo = sqrt(reshape(Phi_Nxz*(price_cum_div(:).^2),[NN Nx Nz])./(price_ex_div.^2) ...
        - exp_ret_stock_t_1mo.^2);
    
    stats.stock_ret_1mo.risk_prem_unlevered = sum(prob_Nxz.*(exp_ret_stock_t_1mo-rf_t_1mo),'all');
    
    var_condE = sum(prob_Nxz.*(exp_ret_stock_t_1mo.^2),'all') - sum(prob_Nxz.*exp_ret_stock_t_1mo,'all')^2;
    E_condvar = sum(prob_Nxz.*(cond_vol_ret_t_1mo.^2),'all');
    stats.stock_ret_1mo.ret_vol_unlevered = sqrt(var_condE + E_condvar);
    
    var_condE = sum(prob_Nxz.*((exp_ret_stock_t_1mo - rf_t_1mo).^2),'all') - sum(prob_Nxz.*(exp_ret_stock_t_1mo - rf_t_1mo),'all')^2;
    stats.stock_ret_1mo.stock_exret_vol_unlevered = sqrt(var_condE + E_condvar);

    if FLAG_RUN_PARALLEL
        
        mean_Xs = zeros(1,4);
        std_Xs = zeros(1,4);
        
        Xs = cell(4,1);
        Xs{1} = log(price_ex_div./gesol.C);
        Xs{2} = log(gesol.C./price_ex_div);
        Xs{3} = gesol.C./price_ex_div;
        Xs{4} = price_ex_div./gesol.C;
        
        parfor i = 1:4
            [mean_Xs(i), std_Xs(i)] = Calc_monthly_vol(prob_Nxz,Xs{i});
        end
        
        stats.stock_ret_1mo.logPDratio_mean = mean_Xs(1);
        stats.stock_ret_1mo.logPDratio_std = std_Xs(1);
        stats.stock_ret_1mo.logDPratio_mean = mean_Xs(2);
        stats.stock_ret_1mo.logDPratio_std = std_Xs(2);
        stats.stock_ret_1mo.DPratio_mean = mean_Xs(3);
        stats.stock_ret_1mo.DPratio_std = std_Xs(3);
        stats.stock_ret_1mo.PDratio_mean = mean_Xs(4);
        stats.stock_ret_1mo.PDratio_std = std_Xs(4);
        
    else
    
        % D = C, P = ex div price
        [stats.stock_ret_1mo.logPDratio_mean,stats.stock_ret_1mo.logPDratio_std] = Calc_monthly_vol(prob_Nxz,log(price_ex_div./gesol.C));
        [stats.stock_ret_1mo.logDPratio_mean,stats.stock_ret_1mo.logDPratio_std] = Calc_monthly_vol(prob_Nxz,log(gesol.C./price_ex_div));
        [stats.stock_ret_1mo.DPratio_mean,stats.stock_ret_1mo.DPratio_std] = Calc_monthly_vol(prob_Nxz,gesol.C./price_ex_div);
        [stats.stock_ret_1mo.PDratio_mean,stats.stock_ret_1mo.PDratio_std] = Calc_monthly_vol(prob_Nxz,price_ex_div./gesol.C);
    
    end
    
    % 1 year returns
    exp_ret_stock_t_1y = ones(NN,Nx,Nz);
    exp_ret_stock_t_1y_sq = ones(NN,Nx,Nz);
    pr = price_ex_div;
    payoff = price_cum_div;
    for n = 1:12
        exp_ret_stock_t_1y(:) = reshape(Phi_Nxz*(payoff(:).*exp_ret_stock_t_1y(:)),[NN Nx Nz])./pr;
        exp_ret_stock_t_1y_sq(:) = reshape(Phi_Nxz*((payoff(:).^2).*exp_ret_stock_t_1y_sq(:)),[NN Nx Nz])./(pr.^2);
    end
    rf_t_1y = 1./gesol.real_bonds.prices(:,:,:,12);
    cond_vol_ret_t_1y = sqrt(exp_ret_stock_t_1y_sq - exp_ret_stock_t_1y.^2);
    
    stats.stock_ret_1y.risk_prem_unlevered = sum(prob_Nxz.*(exp_ret_stock_t_1y-rf_t_1y),'all');
    
    var_condE = sum(prob_Nxz.*(exp_ret_stock_t_1y.^2),'all') - sum(prob_Nxz.*exp_ret_stock_t_1y,'all')^2;
    E_condvar = sum(prob_Nxz.*(cond_vol_ret_t_1y.^2),'all');
    stats.stock_ret_1y.ret_vol_unlevered = sqrt(var_condE + E_condvar);
    
    var_condE = sum(prob_Nxz.*((exp_ret_stock_t_1y - rf_t_1y).^2),'all') - sum(prob_Nxz.*(exp_ret_stock_t_1y - rf_t_1y),'all')^2;
    stats.stock_ret_1y.stock_exret_vol_unlevered = sqrt(var_condE + E_condvar);
    
end